subroutine maxwell_solver
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
      use celeste3d
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!
!     calculate the source for Faraday's law
!     store in jxtildet,jytildet,jztildet
!
!
                  ws = fourpi/clite
                  dtc = clite*cntr*dt
!
                     call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x&
                        , c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y&
                        , c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, &
                        qdnc, jxtildet, jytildet, jztildet)
!
                  call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, &
                     ijkvtx, vol, wate(1,28))
!

                     do n = 1, nvtx
                        ijk = ijkvtx(n)
!
                        rvvol = 1./vvol(ijk)
                        ws = fourpi/clite*vvol(ijk)
!
                        jxtildet(ijk) = (ex0(ijk)-fourpi*jxtildet(ijk)*rvvol*&
                           dtc*dtc+dtc*(jx0(ijk)-ws*jxtilde(ijk))*rvvol)*wate(&
                           ijk,28)
!
                        jytildet(ijk) = (ey0(ijk)-fourpi*jytildet(ijk)*rvvol*&
                           dtc*dtc+dtc*(jy0(ijk)-ws*jytilde(ijk))*rvvol)*wate(&
                           ijk,28)
!
                        jztildet(ijk) = (ez0(ijk)-fourpi*jztildet(ijk)*rvvol*&
                           dtc*dtc+dtc*(jz0(ijk)-ws*jztilde(ijk))*rvvol)*wate(&
                           ijk,28)
!
                     end do
!
                     ex(:ibp2*jbp2*kbp2) = ex0(:ibp2*jbp2*kbp2)
                     ey(:ibp2*jbp2*kbp2) = ey0(:ibp2*jbp2*kbp2)
                     ez(:ibp2*jbp2*kbp2) = ez0(:ibp2*jbp2*kbp2)
!
                     errormaxwell = error
!
                     call maxwell (nvtx, nvtxkm, ijkvtx, iwid, jwid, kwid, &
                        precon, tiny, ncells, ijkcell, itdim, nrg, ibp1, jbp1, &
                        kbp1, wate(1,28), wate(1,29), cdlt, sdlt, strait, dx, dy, dz, &
                        clite, cntr, dt, t, fourpi, itmax, errormaxwell, &
                        periodicx, bx_ini, c1x, &
                        c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, &
                        c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
                        c8z, vvol, vol, x,y,z,qom, qdnv, qdnc, curlx, curly, curlz, wate(1,&
                        1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                        6), wate(1,7), wate(1,8), wate(1,9), wate(1,10), wate(1&
                        ,11), wate(1,12), wate(1,13), wate(1,14), wate(1,15), &
                        wate(1,16), wate(1,17), wate(1,18), wate(1,19), wate(1,&
                        20), wate(1,21), a11, a12, a13, a21, a22, a23, a31, a32&
                        , a33, iphd2, wate(1,22), wate(1,23), wate(1,24), wate(&
                        1,25), wate(1,26), wate(1,27), jxtilde, jytilde, &
                        jztilde, jx0, jy0, jz0, jxtildet, jytildet, jztildet, &
                        pixx, pixy, pixz, piyy, piyz, pizz, divpix, divpiy, &
                        divpiz, colx, coly, colz, udrft, vdrft, wdrft, tsix, &
                        tsiy, tsiz, bxv, byv, bzv, ex0, ey0, ez0, ex, ey, ez, &
                        iter,ex_ext,ey_ext,ez_ext,shock,susxx,susxy&
                        ,susxz,suszx,suszy,suszz)
!
                  citer_maxwell(nh) = iter

         call curlv (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, bxn, byn, bzn, jx0, jy0, jz0)
!

end subroutine maxwell_solver


subroutine maxwell_risolver
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
      use celeste3d
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
real(double) :: qdnv0(itdim,nsp)
!
!     calculate the source for Faraday's law
!     store in jxtildet,jytildet,jztildet
                  ws = fourpi/clite
                  dtc = clite*cntr*dt
!
                     call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x&
                        , c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y&
                        , c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, &
                        qdnc0, jxtildet, jytildet, jztildet)
!
                  call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, &
                     ijkvtx, vol, wate(1,28))
!

                     do n = 1, nvtx
                        ijk = ijkvtx(n)
!
                        rvvol = 1./vvol(ijk)
                        ws = fourpi/clite*vvol(ijk)
!
                        jxtildet(ijk) = (ex0(ijk)-fourpi*jxtildet(ijk)*rvvol*&
                           dtc*dtc+dtc*(jx0(ijk)-ws*j12x(ijk))*rvvol)*wate(&
                           ijk,28)
!
                        jytildet(ijk) = (ey0(ijk)-fourpi*jytildet(ijk)*rvvol*&
                           dtc*dtc+dtc*(jy0(ijk)-ws*j12y(ijk))*rvvol)*wate(&
                           ijk,28)
!
                        jztildet(ijk) = (ez0(ijk)-fourpi*jztildet(ijk)*rvvol*&
                           dtc*dtc+dtc*(jz0(ijk)-ws*j12z(ijk))*rvvol)*wate(&
                           ijk,28)
!
                     end do
!
                     ex(:ibp2*jbp2*kbp2) = ex0(:ibp2*jbp2*kbp2)
                     ey(:ibp2*jbp2*kbp2) = ey0(:ibp2*jbp2*kbp2)
                     ez(:ibp2*jbp2*kbp2) = ez0(:ibp2*jbp2*kbp2)
                     qdnv0(:ibp2*jbp2*kbp2,:nsp) = 0d0
!
                     errormaxwell = error
!
                     call maxwell (nvtx, nvtxkm, ijkvtx, iwid, jwid, kwid, &
                        precon, tiny, ncells, ijkcell, itdim, nrg, ibp1, jbp1, &
                        kbp1, wate(1,28), wate(1,29), cdlt, sdlt, strait, dx, dy, dz, &
                        clite, cntr, dt, t, fourpi, itmax, errormaxwell, &
                        periodicx, bx_ini, c1x, &
                        c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, &
                        c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
                        c8z, vvol, vol, x,y,z,qom, qdnv0, qdnc0, curlx, curly, curlz, wate(1,&
                        1), wate(1,2), wate(1,3), wate(1,4), wate(1,5), wate(1,&
                        6), wate(1,7), wate(1,8), wate(1,9), wate(1,10), wate(1&
                        ,11), wate(1,12), wate(1,13), wate(1,14), wate(1,15), &
                        wate(1,16), wate(1,17), wate(1,18), wate(1,19), wate(1,&
                        20), wate(1,21), a11, a12, a13, a21, a22, a23, a31, a32&
                        , a33, iphd2, wate(1,22), wate(1,23), wate(1,24), wate(&
                        1,25), wate(1,26), wate(1,27), jxtilde, jytilde, &
                        jztilde, jx0, jy0, jz0, jxtildet, jytildet, jztildet, &
                        pixx, pixy, pixz, piyy, piyz, pizz, divpix, divpiy, &
                        divpiz, colx, coly, colz, udrft, vdrft, wdrft, tsix, &
                        tsiy, tsiz, bxv, byv, bzv, ex0, ey0, ez0, ex, ey, ez, &
                        iter,ex_ext,ey_ext,ez_ext,shock,susxx,susxy&
                        ,susxz,suszx,suszy,suszz)
!
                  citer_maxwell(nh) = iter

         call curlv (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, bxn, byn, bzn, jx0, jy0, jz0)
!

end subroutine maxwell_risolver



       subroutine maxwell(nvtx, nvtxkm, ijkvtx, iwid, jwid, kwid, precon, tiny, &
         ncells, ijkcell, itdim, nsp, ibp1, jbp1, kbp1, vvolaxis, divetil, cdlt&
         , sdlt, strait, dx, dy, dz, clite, cntr, dt, t, fourpi, itmax, error &
         , periodicx, bx_ini, &
         c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, &
         c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vvol, vol, x,y,z, qom, qdnv&
         , qdnc, curlx, curly, curlz, exmu, eymu, ezmu, rx, ry, rz, qx, qy, qz, qxm1&
         , qym1, qzm1, aqxm1, aqym1, aqzm1, aqx, aqy, aqz, paqx, paqy, paqz, &
         a11, a12, a13, a21, a22, a23, a31, a32, a33, ijkvtmp, qxtilde, qytilde&
         , qztilde, aqxtilde, aqytilde, aqztilde, jxtilde, jytilde, jztilde, &
         paqxm1, paqym1, paqzm1, gdivx, gdivy, gdivz, qxm2, qym2, qzm2, aqxm2, &
         aqym2, aqzm2, paqxm2, paqym2, paqzm2, qxm3, qym3, qzm3, aqxm3, aqym3, &
         aqzm3, paqxm3, paqym3, paqzm3, bxv, byv, bzv, ex0, ey0, ez0, ex, ey, &
         ez, iter,ex_ext,ey_ext,ez_ext,shock,susxx,susxy,susxz,suszx&
         ,suszy,suszz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nvtx
      integer  :: nvtxkm
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      integer  :: ncells
      integer  :: itdim
      integer  :: nsp
      integer  :: ibp1
      integer  :: jbp1
      integer  :: kbp1
      integer , intent(in) :: itmax
      integer , intent(out) :: iter
      real(double)  :: tiny
      real(double)  :: cdlt
      real(double)  :: sdlt
      real(double)  :: strait
      real(double)  :: dx,dy,dz
      real(double)  :: clite
      real(double)  :: cntr
      real(double)  :: dt
      real(double)  :: t
      real(double)  :: fourpi
      real(double) , intent(in) :: error
      logical , intent(in) :: precon
      logical  :: periodicx
	  real(double) :: ex_ext,ey_ext,ez_ext,bx_ini
      integer  :: ijkvtx(*)
      integer  :: ijkcell(*)
      integer  :: ijkvtmp(*)
      logical :: shock
      real(double)  :: vvolaxis(*)
      real(double)  :: divetil(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double)  :: vvol(*)
      real(double)  :: vol(*)
      real(double)  :: x(*)
      real(double)  :: y(*)
      real(double)  :: z(*)
      real(double)  :: qom(*)
      real(double)  :: qdnv(itdim,*)
      real(double)  :: qdnc(*)
      real(double)  :: curlx(*)
      real(double)  :: curly(*)
      real(double)  :: curlz(*)
      real(double)  :: exmu(*)
      real(double)  :: eymu(*)
      real(double)  :: ezmu(*)
      real(double)  :: rx(*)
      real(double)  :: ry(*)
      real(double)  :: rz(*)
      real(double) , intent(inout) :: qx(*)
      real(double) , intent(inout) :: qy(*)
      real(double) , intent(inout) :: qz(*)
      real(double) , intent(inout) :: qxm1(*)
      real(double) , intent(inout) :: qym1(*)
      real(double) , intent(inout) :: qzm1(*)
      real(double) , intent(inout) :: aqxm1(*)
      real(double) , intent(inout) :: aqym1(*)
      real(double) , intent(inout) :: aqzm1(*)
      real(double)  :: aqx(*)
      real(double)  :: aqy(*)
      real(double)  :: aqz(*)
      real(double)  :: paqx(*)
      real(double)  :: paqy(*)
      real(double)  :: paqz(*)
      real(double)  :: a11(*)
      real(double)  :: a12(*)
      real(double)  :: a13(*)
      real(double)  :: a21(*)
      real(double)  :: a22(*)
      real(double)  :: a23(*)
      real(double)  :: a31(*)
      real(double)  :: a32(*)
      real(double)  :: a33(*)
      real(double)  :: qxtilde(*)
      real(double)  :: qytilde(*)
      real(double)  :: qztilde(*)
      real(double)  :: aqxtilde(*)
      real(double)  :: aqytilde(*)
      real(double)  :: aqztilde(*)
      real(double)  :: jxtilde(*)
      real(double)  :: jytilde(*)
      real(double)  :: jztilde(*)
      real(double) , intent(inout) :: paqxm1(*)
      real(double) , intent(inout) :: paqym1(*)
      real(double) , intent(inout) :: paqzm1(*)
      real(double) , intent(inout) :: gdivx(*)
      real(double) , intent(inout) :: gdivy(*)
      real(double) , intent(inout) :: gdivz(*)
      real(double) , intent(inout) :: qxm2(*)
      real(double) , intent(inout) :: qym2(*)
      real(double) , intent(inout) :: qzm2(*)
      real(double) , intent(inout) :: aqxm2(*)
      real(double) , intent(inout) :: aqym2(*)
      real(double) , intent(inout) :: aqzm2(*)
      real(double) , intent(inout) :: paqxm2(*)
      real(double) , intent(inout) :: paqym2(*)
      real(double) , intent(inout) :: paqzm2(*)
      real(double) , intent(inout) :: qxm3(*)
      real(double) , intent(inout) :: qym3(*)
      real(double) , intent(inout) :: qzm3(*)
      real(double) , intent(inout) :: aqxm3(*)
      real(double) , intent(inout) :: aqym3(*)
      real(double) , intent(inout) :: aqzm3(*)
      real(double) , intent(inout) :: paqxm3(*)
      real(double) , intent(inout) :: paqym3(*)
      real(double) , intent(inout) :: paqzm3(*)
      real(double)  :: bxv(*)
      real(double)  :: byv(*)
      real(double)  :: bzv(*)
      real(double)  :: ex0(*)
      real(double)  :: ey0(*)
      real(double)  :: ez0(*)
      real(double)  :: ex(*)
      real(double)  :: ey(*)
      real(double)  :: ez(*)
      real(double)  :: susxx(*)
      real(double)  :: susxy(*)
      real(double)  :: susxz(*)
      real(double)  :: suszx(*)
      real(double)  :: suszy(*)
      real(double)  :: suszz(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ibp2, jbp2, kbp2, n, ijk, k, j, i, k2, ijk2, i2, nvtmp
      integer :: istart, jstart, kstart,ngmres
      integer :: ijkgmres(itdim)
      real(double) :: jx0, jy0, jz0, lambda, lambm1, jxtildet, jytildet, &
         jztildet, lambm2, lambm3, dnorm, rsum, bottom, botm1, botm2, botm3, &
         dero, dern, derd, wsqx, wsqy, wsqz, wsaqx, wsaqy, wsaqz, rnorm, alfa, &
         dxnorm, xnorm,factor, norm_in, norm_b
      logical :: testgl, precongl
!      real(double):: susxx(ijkvtx(nvtx)), susxy(ijkvtx(nvtx)), susxz(ijkvtx(nvtx)),&
!           suszx(ijkvtx(nvtx)), suszy(ijkvtx(nvtx))&
!           ,suszz(ijkvtx(nvtx)),susyx(ijkvtx(nvtx)),susyy(ijkvtx(nvtx)),susyz(ijkvtx(nvtx))&
!-----------------------------------------------
!
!
!
!
      ibp2 = ibp1 + 1
      jbp2 = jbp1 + 1
      kbp2 = kbp1 + 1
!
      istart=3
      jstart=3
      kstart=3
!
      if (periodicx) istart=2
      if (periodicy) jstart=2
      if (periodicz) kstart=2
!
      call list (istart, ibp1, jstart, jbp1, kstart, kbp1, iwid, jwid, kwid, ngmres, ijkgmres)
!	write(*,*)'number of unknowns',ngmres
!
      call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, ngmres, ijkgmres, vol, &
         vvolaxis)
!
!     calculate A*E0
!
!     ****************************************************************
!
!     CONSTRUCT JACOBI PRECONDITIONER
!
!     ***************************************************************
!
!
      if (precon) then
!
         call block (ngmres, ijkgmres, nsp, itdim, iwid, jwid, kwid, c1x, c2x, c3x&
            , c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, &
            c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vvolaxis, vol, vvol, clite&
            , dt, qom, qdnv, cntr, bxv, byv, bzv, exmu, eymu, ezmu, a11, a12, &
            a13, a21, a22, a23, a31, a32, a33)
!
!
         call ludecomp (ngmres, ijkgmres, a11, a12, a13, a21, a22, a23, a31, a32, &
            a33)
!
      else
!
         dnorm = 1.d0
!dir$ ivdep
         do n = 1, ngmres
!
            ijk = ijkgmres(n)
!
            a11(ijk) = dnorm
            a12(ijk) = 0.d0
            a13(ijk) = 0.d0
!cc
            a21(ijk) = 0.d0
            a22(ijk) = dnorm
            a23(ijk) = 0.d0
!ccc
            a31(ijk) = 0.d0
            a32(ijk) = 0.d0
            a33(ijk) = dnorm
!
         end do
!
      endif
!
!     END PRE-CONDITIIONER
!
!     ********************************************************************
!
!     CALCULATE INITIAL RESIDUE
!
!     **********************************************************************
!
	factor=1d0
!
        call sustensor(nvtx, ijkvtx, nsp, itdim, qdnv, qom, bxv, byv, bzv, dt&
         , 1.d0, clite, cntr, susxx, susxy, susxz,suszx, suszy,suszz)
!
         call   bc_maxwell(t,ibp2,jbp2,kbp2,nsp, dx,dy,dz, dt,ex0, ey0, ez0,  &
         bxv, byv, bzv, ex, ey, ez, vvol,x,y,z,periodicx,ex_ext&
         ,ey_ext,ez_ext,bx_ini,shock,nvtx,ijkvtx,itdim,qdnv,qom,&
         clite,cntr,jxtilde,jztilde,susxx,susxy,susxz,suszx,suszy&
         ,suszz)
         call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
            cdlt, ex, ey, ez, periodicx)
!
      call avec (nvtmp, nvtxkm, ijkvtmp, iwid, jwid, kwid, tiny, ncells, ijkcell&
         , qom, qdnv, nsp, itdim, vvol, ibp1, jbp1, kbp1, vvolaxis, cdlt, sdlt&
         , strait, dx, dy, dz, periodicx, periodicz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, &
         c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
         c8z, vol, curlx, curly, curlz, exmu, eymu, ezmu, bxv, byv, bzv, cntr, &
         dt, 1.d0, clite, fourpi, ex, ey, ez, rx, ry, rz,factor)

!dir$ ivdep
      do n = 1, ngmres
!
         ijk = ijkgmres(n)
!
!
         rx(ijk) = gdivx(ijk) - rx(ijk)
         ry(ijk) = gdivy(ijk) - ry(ijk)
         rz(ijk) = gdivz(ijk) - rz(ijk)
!
      end do
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         rx, ry, rz, periodicx)

!
!
!     END CALCULATION OF RESIDUE
!
!     ****************************************************************
!
      qxtilde(:ibp2*jbp2*kbp2) = ex(:ibp2*jbp2*kbp2)
      qytilde(:ibp2*jbp2*kbp2) = ey(:ibp2*jbp2*kbp2)
      qztilde(:ibp2*jbp2*kbp2) = ez(:ibp2*jbp2*kbp2)
!
!
      rsum = 0.0d0
! ijkgmx=0
!dir$ ivdep
      do n = 1, ngmres
!
         ijk = ijkgmres(n)
!
         qx(ijk) = 0.0d0
         aqx(ijk) = 0.0d0
         paqx(ijk) = 0.0d0
!
         qy(ijk) = 0.0d0
         aqy(ijk) = 0.0d0
         paqy(ijk) = 0.0d0
!
         qz(ijk) = 0.0d0
         aqz(ijk) = 0.0d0
         paqz(ijk) = 0.0d0
!
         qxm1(ijk) = 0.0d0
         aqxm1(ijk) = 0.0d0
         paqxm1(ijk) = 0.0d0
!
         qym1(ijk) = 0.0d0
         aqym1(ijk) = 0.0d0
         paqym1(ijk) = 0.0d0
!
         qzm1(ijk) = 0.0d0
         aqzm1(ijk) = 0.0d0
         paqzm1(ijk) = 0.0d0
!
         qxm2(ijk) = 0.0d0
         aqxm2(ijk) = 0.0d0
         paqxm2(ijk) = 0.0d0
!
         qym2(ijk) = 0.0d0
         aqym2(ijk) = 0.0d0
         paqym2(ijk) = 0.0d0
!
         qzm2(ijk) = 0.0d0
         aqzm2(ijk) = 0.0d0
         paqzm2(ijk) = 0.0d0
!
         qxm3(ijk) = 0.0d0
         aqxm3(ijk) = 0.0d0
         paqxm3(ijk) = 0.0d0
!
         qym3(ijk) = 0.0d0
         aqym3(ijk) = 0.0d0
         paqym3(ijk) = 0.0d0
!
         qzm3(ijk) = 0.0d0
         aqzm3(ijk) = 0.0d0
         paqzm3(ijk) = 0.0d0
!
!
         rsum = rsum + abs(gdivx(ijk)) + abs(gdivy(ijk)) + abs(gdivz(ijk))
!
!
!
      end do
!
      if (rsum < 1d-10) then
         iter = 0
         go to 10
      endif
!
      bottom = 0.0d0
      botm1 = 0.0d0
      botm2 = 0.0d0
      botm3 = 0.0d0
!
!
!     **************  MAIN ITERATION LOOP  **************************
!     begin conjugate gradient interation
!
      do iter = 1, itmax
!
!
!     APPLY BLOCK PRE-CONDITIONING
!
!
         call precond (ngmres, ijkgmres, a11, a12, a13, a21, a22, a23, a31, a32, &
            a33, rx, ry, rz, qxtilde, qytilde, qztilde)
!
!
!
!dir$ ivdep
         do n = 1, ngmres
!
            ijk = ijkgmres(n)
!
            qxtilde(ijk) = ex(ijk) + qxtilde(ijk)
            qytilde(ijk) = ey(ijk) + qytilde(ijk)
            qztilde(ijk) = ez(ijk) + qztilde(ijk)
!
         end do
!
        call sustensor(nvtx, ijkvtx, nsp, itdim, qdnv, qom, bxv, byv, bzv, dt&
         , 1.d0, clite, cntr, susxx, susxy, susxz,suszx, suszy,suszz)
!
         call   bc_maxwell(t,ibp2,jbp2,kbp2,nsp, dx,dy,dz, dt,ex0, ey0, ez0,  &
         bxv, byv, bzv, qxtilde, qytilde, qztilde, vvol,x,y,z&
         ,periodicx,ex_ext,ey_ext,ez_ext,bx_ini,shock,nvtx,ijkvtx,itdim,qdnv,qom,&
         clite,cntr,jxtilde,jztilde,susxx,susxy,susxz,suszx,suszy,suszz)
         call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
            cdlt, qxtilde, qytilde, qztilde, periodicx)
!
         call avec (nvtmp, nvtxkm, ijkvtmp, iwid, jwid, kwid, tiny, ncells, &
            ijkcell, qom, qdnv, nsp, itdim, vvol, ibp1, jbp1, kbp1, vvolaxis, &
            cdlt, sdlt, strait, dx, dy, dz, periodicx, periodicz, c1x, c2x, c3x, c4x, c5x, c6x, &
            c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
            c4z, c5z, c6z, c7z, c8z, vol, curlx, curly, curlz, exmu, eymu, ezmu&
            , bxv, byv, bzv, cntr, dt, 1., clite, fourpi, qxtilde, qytilde, &
            qztilde, aqxtilde, aqytilde, aqztilde,factor)

            do n = 1, ngmres
!
               ijk = ijkgmres(n)
!
               aqxtilde(ijk) = rx(ijk) - (gdivx(ijk)-aqxtilde(ijk))
               aqytilde(ijk) = ry(ijk) - (gdivy(ijk)-aqytilde(ijk))
               aqztilde(ijk) = rz(ijk) - (gdivz(ijk)-aqztilde(ijk))
!
               qxtilde(ijk) = qxtilde(ijk) - ex(ijk)
               qytilde(ijk) = qytilde(ijk) - ey(ijk)
               qztilde(ijk) = qztilde(ijk) - ez(ijk)
!
            end do
!
         lambda = 0.0d0
         lambm1 = 0.0d0
         lambm2 = 0.0d0
         lambm3 = 0.0d0
!
!dir$ ivdep
         do n = 1, ngmres
!
            ijk = ijkgmres(n)
!
            lambda = lambda + aqxtilde(ijk)*paqx(ijk) + aqytilde(ijk)*paqy(ijk)&
                + aqztilde(ijk)*paqz(ijk)
!
            lambm1 = lambm1 + aqxtilde(ijk)*paqxm1(ijk) + aqytilde(ijk)*paqym1(&
               ijk) + aqztilde(ijk)*paqzm1(ijk)
!
            lambm2 = lambm2 + aqxtilde(ijk)*paqxm2(ijk) + aqytilde(ijk)*paqym2(&
               ijk) + aqztilde(ijk)*paqzm2(ijk)
!
            lambm3 = lambm3 + aqxtilde(ijk)*paqxm3(ijk) + aqytilde(ijk)*paqym3(&
               ijk) + aqztilde(ijk)*paqzm3(ijk)
!
!
         end do
!
         lambda = lambda/(bottom + 1.d-30)
         lambm1 = lambm1/(botm1 + 1.d-30)
         lambm2 = lambm2/(botm2 + 1.d-30)
         lambm3 = lambm3/(botm3 + 1.d-30)
!
!dir$ ivdep
         do n = 1, ngmres
!
            ijk = ijkgmres(n)
!
!
            wsqx = qxtilde(ijk) - lambda*qx(ijk) - lambm1*qxm1(ijk) - lambm2*&
               qxm2(ijk) - lambm3*qxm3(ijk)
            wsqy = qytilde(ijk) - lambda*qy(ijk) - lambm1*qym1(ijk) - lambm2*&
               qym2(ijk) - lambm3*qym3(ijk)
            wsqz = qztilde(ijk) - lambda*qz(ijk) - lambm1*qzm1(ijk) - lambm2*&
               qzm2(ijk) - lambm3*qzm3(ijk)
!
            wsaqx = aqxtilde(ijk) - lambda*aqx(ijk) - lambm1*aqxm1(ijk) - &
               lambm2*aqxm2(ijk) - lambm3*aqxm3(ijk)
            wsaqy = aqytilde(ijk) - lambda*aqy(ijk) - lambm1*aqym1(ijk) - &
               lambm2*aqym2(ijk) - lambm3*aqym3(ijk)
            wsaqz = aqztilde(ijk) - lambda*aqz(ijk) - lambm1*aqzm1(ijk) - &
               lambm2*aqzm2(ijk) - lambm3*aqzm3(ijk)
!
            qxm3(ijk) = qxm2(ijk)
            aqxm3(ijk) = aqxm2(ijk)
            paqxm3(ijk) = paqxm2(ijk)
!
            qym3(ijk) = qym2(ijk)
            aqym3(ijk) = aqym2(ijk)
            paqym3(ijk) = paqym2(ijk)
!
            qzm3(ijk) = qzm2(ijk)
            aqzm3(ijk) = aqzm2(ijk)
            paqzm3(ijk) = paqzm2(ijk)
!
            qxm2(ijk) = qxm1(ijk)
            aqxm2(ijk) = aqxm1(ijk)
            paqxm2(ijk) = paqxm1(ijk)
!
            qym2(ijk) = qym1(ijk)
            aqym2(ijk) = aqym1(ijk)
            paqym2(ijk) = paqym1(ijk)
!
            qzm2(ijk) = qzm1(ijk)
            aqzm2(ijk) = aqzm1(ijk)
            paqzm2(ijk) = paqzm1(ijk)
!
            qxm1(ijk) = qx(ijk)
            aqxm1(ijk) = aqx(ijk)
            paqxm1(ijk) = paqx(ijk)
!
            qym1(ijk) = qy(ijk)
            aqym1(ijk) = aqy(ijk)
            paqym1(ijk) = paqy(ijk)
!
            qzm1(ijk) = qz(ijk)
            aqzm1(ijk) = aqz(ijk)
            paqzm1(ijk) = paqz(ijk)
!
            qx(ijk) = wsqx
            aqx(ijk) = wsaqx
!
            qy(ijk) = wsqy
            aqy(ijk) = wsaqy
!
            qz(ijk) = wsqz
            aqz(ijk) = wsaqz
!
         end do
!
         call precond (ngmres, ijkgmres, a11, a12, a13, a21, a22, a23, a31, a32, &
            a33, aqx, aqy, aqz, paqx, paqy, paqz)
!
!
!
         rnorm = 0.0d0
         alfa = 0.0d0
         bottom = 0.0d0
         botm1 = 0.0d0
         botm2 = 0.0d0
         botm3 = 0.0d0
!
!
!
!dir$ ivdep
         do n = 1, ngmres
!
            ijk = ijkgmres(n)
!
!
!           rnorm = rnorm + abs(rx(ijk)) + abs(ry(ijk)) + abs(rz(ijk))
!
            alfa = alfa + rx(ijk)*paqx(ijk) + ry(ijk)*paqy(ijk) + rz(ijk)*paqz(&
               ijk)
!
            bottom = bottom + aqx(ijk)*paqx(ijk) + aqy(ijk)*paqy(ijk) + aqz(ijk&
               )*paqz(ijk)
!
            botm1 = botm1 + aqxm1(ijk)*paqxm1(ijk) + aqym1(ijk)*paqym1(ijk) + &
               aqzm1(ijk)*paqzm1(ijk)
!
            botm2 = botm2 + aqxm2(ijk)*paqxm2(ijk) + aqym2(ijk)*paqym2(ijk) + &
               aqzm2(ijk)*paqzm2(ijk)
!
            botm3 = botm3 + aqxm3(ijk)*paqxm3(ijk) + aqym3(ijk)*paqym3(ijk) + &
               aqzm3(ijk)*paqzm3(ijk)
!
!
         end do
!
         alfa = alfa/(bottom + 1.d-30)
!
!dir$ ivdep
            do n = 1, ngmres
!
               ijk = ijkgmres(n)
!
               ex(ijk) = ex(ijk) + alfa*qx(ijk)
               ey(ijk) = ey(ijk) + alfa*qy(ijk)
               ez(ijk) = ez(ijk) + alfa*qz(ijk)
!
               rx(ijk) = rx(ijk) - alfa*aqx(ijk)
               ry(ijk) = ry(ijk) - alfa*aqy(ijk)
               rz(ijk) = rz(ijk) - alfa*aqz(ijk)
!
            end do
!
         dxnorm = 0.0
         xnorm = 0.0
!
!dir$ ivdep
            do n = 1, ngmres
!
               ijk = ijkgmres(n)
!
               rnorm = rnorm + abs(rx(ijk)) + abs(ry(ijk)) + abs(rz(ijk))
!
               dxnorm = dxnorm + abs(qx(ijk)) + abs(qy(ijk)) + abs(qz(ijk))
!
               xnorm = xnorm + abs(ex(ijk)) + abs(ey(ijk)) + abs(ez(ijk))
!
            end do
!
         dxnorm = dxnorm*abs(alfa)
!
!      write(*,*)'Iter',iter
!      write(*,*)'TEST Conv2','dxnorm=',dxnorm,'xnorm=',xnorm
!      write(*,*)'TEST Conv2','rnorm=',rnorm,'rsum=',rsum
         if (dxnorm<=error*xnorm .and. rnorm<=error*rsum) go to 10
!
         call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, &
            cdlt, rx, ry, rz, periodicx)
!
!
!
      end do
!
      write (*, *) 'maxwell did not converge', ngmres,rnorm,rsum,dxnorm,xnorm
!
         ex(:itdim)=0d0
         ey(:itdim)=0d0
         ez(:itdim)=0d0
!
        call sustensor(nvtx, ijkvtx, nsp, itdim, qdnv, qom, bxv, byv, bzv, dt&
         , 1.d0, clite, cntr, susxx, susxy, susxz,suszx, suszy,suszz)
!
      call   bc_maxwell(t,ibp2,jbp2,kbp2,nsp, dx,dy,dz, dt,ex0, ey0, ez0,  &
         bxv, byv, bzv, ex, ey, ez, vvol,x,y,z,periodicx,ex_ext&
         ,ey_ext,ez_ext,bx_ini,shock,nvtx,ijkvtx,itdim,qdnv,qom,&
         clite,cntr,jxtilde,jztilde,susxx,susxy,susxz,suszx,suszy&
         ,suszz)

      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         ex, ey, ez, periodicx)
!
!          write(*,*)'failure',sum(abs(rx(ijkgmres(1:ngmres)))) ,sum(abs(ry(ijkgmres(1:ngmres)))),sum(abs(rz(ijkgmres(1:ngmres))))
!         do i=2,ibp1
!         do j=2,jbp1
!         do k=2,kbp1
!         ijk=1+iwid*(i-1)+jwid*(j-1)+kwid*(k-1)
!         write(666,*) gdivy(ijk),ey(ijk),ry(ijk)
!         enddo
!         enddo
!         enddo
!
!
      return
!     **********************End of Iteration **************************
!
   10 continue
      write (*, *) 'maxwell converged', iter, rnorm
!
!         ex(:itdim)=gdivx(:itdim)
!         ey(:itdim)=gdivy(:itdim)
!         ez(:itdim)=gdivz(:itdim)
        call sustensor(nvtx, ijkvtx, nsp, itdim, qdnv, qom, bxv, byv, bzv, dt&
         , 1.d0, clite, cntr, susxx, susxy, susxz,suszx, suszy,suszz)
!
      call   bc_maxwell(t,ibp2,jbp2,kbp2,nsp, dx,dy,dz, dt,ex0, ey0, ez0,   &
         bxv, byv, bzv, ex, ey, ez, vvol,x,y,z,periodicx,ex_ext&
         ,ey_ext,ez_ext,bx_ini,shock,nvtx,ijkvtx,itdim,qdnv,qom,&
         clite,cntr,jxtilde,jztilde,susxx,susxy,susxz,suszx,suszy&
         ,suszz)

      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         ex, ey, ez, periodicx)
!
!
!     ********************** END MAXWELL **********************
!
      return
      end subroutine maxwell


      subroutine avec(nvtmp, nvtxkm, ijkvtmp, iwid, jwid, kwid, tiny, ncells, &
         ijkcell, qom, qdnv, nsp, itdim, vvol, ibp1, jbp1, kbp1, vvolaxis, cdlt&
         , sdlt, strait, dx,dy,dz, periodicx, periodicz,c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x&
         , c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z&
         , c7z, c8z, vol, curlvx, curlvy, curlvz, exmu, eymu, ezmu, bxv, byv, &
         bzv, cntr, dt, transpos, clite, fourpi, ex, ey, ez, hx, hy, hz,factor)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nvtmp
      integer  :: nvtxkm
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      integer  :: ncells
      integer  :: nsp
      integer  :: itdim
      integer  :: ibp1
      integer  :: jbp1
      integer  :: kbp1
      real(double)  :: tiny
      real(double)  :: factor
      real(double)  :: cdlt
      real(double)  :: sdlt
      real(double)  :: strait
      real(double)  :: dx,dy,dz
      real(double) , intent(in) :: cntr
      real(double)  :: dt
      real(double)  :: transpos
      real(double)  :: clite
      real(double)  :: fourpi
      logical  :: periodicx, periodicz
      integer  :: ijkvtmp(*)
      integer  :: ijkcell(*)
      real(double)  :: qom(*)
      real(double)  :: qdnv(itdim,*)
      real(double) , intent(in) :: vvol(*)
      real(double) , intent(in) :: vvolaxis(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double)  :: vol(*)
      real(double)  :: curlvx(*)
      real(double)  :: curlvy(*)
      real(double)  :: curlvz(*)
      real(double)  :: exmu(*)
      real(double)  :: eymu(*)
      real(double)  :: ezmu(*)
      real(double)  :: bxv(*)
      real(double)  :: byv(*)
      real(double)  :: bzv(*)
      real(double)  :: ex(*)
      real(double)  :: ey(*)
      real(double)  :: ez(*)
      real(double) , intent(out) :: hx(*)
      real(double) , intent(out) :: hy(*)
      real(double) , intent(out) :: hz(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j, i, ijk, n, ncellstmp
	  integer, dimension(itdim) :: ijkcelltmp
      real(double) :: rvvol, dtsq, dtc, dtcsq
	  real(double), dimension(itdim) :: divD
	  logical :: cartesian
!-----------------------------------------------
!
      cartesian = .TRUE.
!
!     calculate the curl of E at cell centers and store in curlc
!
!     AVEC calls CULRC, CURLV, and MUDOTE
!
      call list (1, ibp1, 1, jbp1, 1, kbp1, iwid, jwid, kwid, nvtmp, ijkvtmp)
!
!
!     EX COMPONENT
!
      call gradc (nvtmp, ijkvtmp, kbp1, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, vol, ex, exmu, eymu, ezmu)
!
!     calculate the divergence of gradc at cell vertices and store in curlv
!
      call divv (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, exmu, eymu, ezmu, curlvx)
!
!     ey component
!
      call gradc (nvtmp, ijkvtmp, kbp1, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, vol, ey, exmu, eymu, ezmu)
!
!     calculate the divergence of gradc at cell vertices and store in curlv
!
      call divv (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, exmu, eymu, ezmu, curlvy)
!
!
!     EZ COMPONENT
!
      call gradc (nvtmp, ijkvtmp, kbp1, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, vol, ez, exmu, eymu, ezmu)
!
!     calculate the divergence of gradc at cell vertices and store in curlv
!
      call divv (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, exmu, eymu, ezmu, curlvz)
!
      call list (2, ibp1 + 1, 2, jbp1 + 1, 2, kbp1 + 1, iwid, jwid, kwid, nvtmp&
         , ijkvtmp)
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         curlvx, curlvy, curlvz, periodicx)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, curlvx, curlvy, curlvz)
!
!     normalize axis divergence
!
      if (.not.cartesian) then
         do j = 2, jbp1 + 1
            do i = 2, ibp1 + 1
!
               ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + kwid
               rvvol = 1./vvol(ijk)
               curlvx(ijk) = curlvx(ijk)*vvolaxis(ijk)*rvvol
               curlvy(ijk) = curlvy(ijk)*vvolaxis(ijk)*rvvol
               curlvz(ijk) = curlvz(ijk)*vvolaxis(ijk)*rvvol
!
!
            end do
         end do
      endif
!
!
!
!     project susceptibility tensor on to E
!
      call mudote (nvtmp, ijkvtmp, nsp, itdim, qdnv, qom, bxv, byv, bzv, dt, &
         transpos, clite, ex, ey, ez, exmu, eymu, ezmu)
!      write(*,*)'exmu',sum(exmu(ijkvtmp(:nvtmp))**2),sum(eymu(ijkvtmp(:nvtmp))**2),sum(ezmu(ijkvtmp(:nvtmp))**2)
!
      call graddiv (ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1, kbp1, &
         cdlt, sdlt, periodicx, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, &
         c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
         c8z, vol, exmu, eymu, ezmu, divD, hx, hy, hz)
!
      dtsq = dt**2
      dtc = clite*cntr*dt
      dtcsq = dtc**2
!
!	hx(:itdim)=0d0
!	hy(:itdim)=0d0
!	hz(:itdim)=0d0
!
     if (cartesian) then
!dir$ ivdep
         do n = 1, nvtmp
!

            ijk = ijkvtmp(n)
!
            hx(ijk) = (-dtcsq*(curlvx(ijk)+0.5*cntr*dtsq*hx(ijk))) + vvol(ijk)*&
               (0.5*cntr*dtsq*exmu(ijk)+ex(ijk))
            hy(ijk) = (-dtcsq*(curlvy(ijk)+0.5*cntr*dtsq*hy(ijk))) + vvol(ijk)*&
               (0.5*cntr*dtsq*eymu(ijk)+ey(ijk))
            hz(ijk) = (-dtcsq*(curlvz(ijk)+0.5*cntr*dtsq*hz(ijk))) + vvol(ijk)*&
               (0.5*cntr*dtsq*ezmu(ijk)+ez(ijk))
!
         end do
!
      else
!dir$ ivdep
         do n = 1, nvtmp
!

            ijk = ijkvtmp(n)
!
            hx(ijk) = (-dtcsq*curlvx(ijk)) + vvolaxis(ijk)*(0.5*cntr*dtsq*exmu(&
               ijk)+ex(ijk))
            hy(ijk) = (-dtcsq*curlvy(ijk)) + vvolaxis(ijk)*(0.5*cntr*dtsq*eymu(&
               ijk)+ey(ijk))
            hz(ijk) = (-dtcsq*curlvz(ijk)) + vvolaxis(ijk)*(0.5*cntr*dtsq*ezmu(&
               ijk)+ez(ijk))
!
         end do
      endif
!
!      	hx(:itdim)=ex(:itdim)
!	hy(:itdim)=ey(:itdim)
!	hz(:itdim)=ez(:itdim)
!
      return
      end subroutine avec


      subroutine axisgrad(ibp1, jbp1, iwid, jwid, kwid, gradxf, gradyf, gradzf)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ibp1
      integer  :: jbp1
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      real(double)  :: gradxf(*)
      real(double)  :: gradyf(*)
      real(double)  :: gradzf(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j, i, ijk
      real(double) :: gradx, grady, gradz
!-----------------------------------------------
!
!
!     calculate gradient at the axis by averaging gradients in the poloidal angle
!
!     WARNING!!! THIS VERSION ONLY WORKS FOR CARTESIAN GEOETRIES
      return
!
      end subroutine axisgrad




      subroutine axisvol(ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, ijkvtx, vol&
         , vvolaxis)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ibp1
      integer , intent(in) :: jbp1
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      integer , intent(in) :: nvtx
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(in) :: vvol(*)
      real(double) , intent(in) :: vol(*)
      real(double) , intent(out) :: vvolaxis(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, j, i, ijk
!-----------------------------------------------
!
      vvolaxis(ijkvtx(:nvtx)) = vvol(ijkvtx(:nvtx))
!
!
!     WARNING!!! THIS VERSION ONLY WORKS FOR CARTESIAN GEOETRIES
return
!
      do j = 2, jbp1 + 1
!
         vvolaxis(iwid+1+(j-1)*jwid+kwid:ibp1*iwid+1+(j-1)*jwid+kwid:iwid) = &
            0.125*(vol(iwid+1+(j-1)*jwid+kwid:ibp1*iwid+1+(j-1)*jwid+kwid:iwid)&
            +vol(1+(j-1)*jwid+kwid:(ibp1-1)*iwid+1+(j-1)*jwid+kwid:iwid)+vol(1+&
            (j-2)*jwid+kwid:(ibp1-1)*iwid+1+(j-2)*jwid+kwid:iwid)+vol(1+iwid+(j&
            -2)*jwid+kwid:ibp1*iwid+1+(j-2)*jwid+kwid:iwid))
!
      end do
!
      return
      end subroutine axisvol

      subroutine bc_maxwell(t,nxp,nyp,nzp,nsp, dx,dy,dz, dt,ex0, ey0, ez0,   &
         bx, by, bz, ex, ey, ez, vvol,x,y,z,periodicx,ex_ext,ey_ext&
         ,ez_ext,bx_ini,shock, nvtx,ijkvtx,itdim,qdnv,qom,&
         clite,cntr,jxtilde,jztilde,susxx,susxy,susxz,suszx,suszy&
         ,suszz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary

! Master method to control application of boundary conditions at each face of the system in 3D
! The conditions are applied to all first and last active node in each direction
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nxp,nyp,nzp
      integer  :: nsp
      logical  :: periodicx,shock
      real(double)  :: dx,dy,dz,dt,t,bx_ini,ex_ext,ey_ext,ez_ext
      real(double)  :: qom(nsp)
      real(double)  :: ex0(nxp, nyp, nzp)
      real(double)  :: ey0(nxp, nyp, nzp)
      real(double)  :: ez0(nxp, nyp, nzp)
      real(double)  :: bx(nxp, nyp, nzp)
      real(double)  :: by(nxp, nyp, nzp)
      real(double)  :: bz(nxp, nyp, nzp)
      real(double)  :: ex(nxp, nyp, nzp)
      real(double)  :: ey(nxp, nyp, nzp)
      real(double)  :: ez(nxp, nyp, nzp)
      real(double)  :: vvol(nxp, nyp, nzp)
      real(double)  :: x(nxp, nyp, nzp)
      real(double)  :: y(nxp, nyp, nzp)
      real(double)  :: z(nxp, nyp, nzp)
      real(double)  :: susxx(nxp, nyp, nzp)
      real(double)  :: susxy(nxp, nyp, nzp)
      real(double)  :: susxz(nxp, nyp, nzp)
      real(double)  :: suszx(nxp, nyp, nzp)
      real(double)  :: suszy(nxp, nyp, nzp)
      real(double)  :: suszz(nxp, nyp, nzp)
      real(double)  :: jxtilde(nxp, nyp, nzp)
      real(double)  :: jztilde(nxp, nyp, nzp)
      real(double), intent(in)  :: qdnv,clite,cntr
      integer, intent(in) :: nvtx, itdim
      integer, intent(in) :: ijkvtx(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: i,j,k,ibar
      real(double) :: a_m,o_m,g_x,v_int,vavg,bavg2,timescale,z_c,pi
!-----------------------------------------------
!-----------------------------------------------
!
! ... to find the elements of the susceptibility tensor
!      write(*,*)'in ... susxx(2,2,2) is  ', susxx(2,2,2)

!
timescale=bx_ini
! the standard for the Newton challenge is a_m=2
!a_m=2.d0*1.5d0
ibar=nxp-2
a_m=ey_ext
o_m=.05d0*timescale
pi=acos(-1.d0)

!   Boundaries in X
!
!   min boundary in x
!

	i=2
      do j=2,nyp
	  do k=2,nzp
	  if (periodicx) then

	  ! do nothing

	  else

	  if(bcMx.eq.1) then
	    ! apply conductor
        ey(i,j,k)=0d0
        ez(i,j,k)=0d0
        ex(i,j,k) = (ex0(i,j,k)-susxy(i,j,k)*ey(i,j,k)-susxz(i,j,k)&
                *ez(i,j,k) -jxtilde(i,j,k)*cntr*dt*1.d0)/susxx(i,j,k)*1.d0
	  elseif(bcMx.eq.2) then
	    !apply mirror ... and open
	    ex(i,j,k)=0d0
	    ey(i,j,k)=(4d0*ey(i+1,j,k)-ey(i+2,j,k))/3d0
	    ez(i,j,k)=(4d0*ez(i+1,j,k)-ez(i+2,j,k))/3d0
      elseif(bcMx.eq.0) then
        !apply open, OURS
        ey(i,j,k)=(4d0*ey(i+1,j,k)-ey(i+2,j,k))/3d0 ! if Ey continous, Bz would always be 0 at boundary
        ez(i,j,k)=(4d0*ez(i+1,j,k)-ez(i+2,j,k))/3d0  ! Ez continuous
        ex(i,j,k) = (ex0(i,j,k)-susxy(i,j,k)*ey(i,j,k)-susxz(i,j,k)&
                *ez(i,j,k) -jxtilde(i,j,k)*cntr*dt*1.d0)/susxx(i,j,k)*0.d0 ! Ex from local density
      elseif(bcMx.eq.-1) then
        !apply open, Divin07
        ey(i,j,k)=(4d0*ey(i+1,j,k)-ey(i+2,j,k))/3d0
        ez(i,j,k)=0d0
        ex(i,j,k)=(4d0*ex(i+1,j,k)-ex(i+2,j,k))/3d0
      elseif(bcMx.eq.-2) then
        !apply open, Horiuchi01
        ey(i,j,k)=(5d0*ey(i+1,j,k)-4d0*ey(i+2,j,k)+ey(i+3,j,k))/2d0 !dEy/dx continous
        ez(i,j,k)=(4d0*ez(i+1,j,k)-ez(i+2,j,k))/3d0  ! Ez continuous
        ex(i,j,k)=(4d0*ex(i+1,j,k)-ex(i+2,j,k))/3d0  ! Ex continuous
      else
        !apply open ..the Pritchett's version .. actually it's not a real open BC
        ey(i,j,k)=(4d0*ey(i+1,j,k)-ey(i+2,j,k))/3d0
        ez(i,j,k)=0d0
        ex(i,j,k)= (ex0(i,j,k)-susxy(i,j,k)*ey(i,j,k)-susxz(i,j,k)&
               *ez(i,j,k) -jxtilde(i,j,k)*cntr*dt*1.d0)/susxx(i,j,k)*1.d0
	  end if

	  end if
	  enddo
	  enddo
!
!     max boundary in x
!
	i=nxp
      do j=2,nyp
	  do k=2,nzp
	  if (periodicx) then
	  ex(i,j,k)=ex(2,j,k)
  	  ey(i,j,k)=ey(2,j,k)
  	  ez(i,j,k)=ez(2,j,k)
	  else
	  if(bcMx.eq.1) then
	    ! apply conductor
        ey(i,j,k)=0d0
        ez(i,j,k)=0d0
        ex(i,j,k) = (ex0(i,j,k)-susxy(i,j,k)*ey(i,j,k)-susxz(i,j,k)&
               *ez(i,j,k) -jxtilde(i,j,k)*cntr*dt*1.d0)/susxx(i,j,k)*1.d0
	  elseif(bcMx.eq.2) then
	    !apply mirror ... and open
	    ex(i,j,k)=0d0
	    ey(i,j,k)=(4d0*ey(i-1,j,k)-ey(i-2,j,k))/3d0
	    ez(i,j,k)=(4d0*ez(i-1,j,k)-ez(i-2,j,k))/3d0
      elseif(bcMx.eq.0) then
        !apply open, OURS
        ey(i,j,k)=(4d0*ey(i-1,j,k)-ey(i-2,j,k))/3d0
        ez(i,j,k)=(4d0*ez(i-1,j,k)-ez(i-2,j,k))/3d0
        ex(i,j,k) = (ex0(i,j,k)-susxy(i,j,k)*ey(i,j,k)-susxz(i,j,k)&
               *ez(i,j,k) -jxtilde(i,j,k)*cntr*dt*1.d0)/susxx(i,j,k)*0.d0
      elseif(bcMx.eq.-1) then
        !apply open, Divin07
        ey(i,j,k)=(4d0*ey(i-1,j,k)-ey(i-2,j,k))/3d0
        ez(i,j,k)=0d0
        ex(i,j,k)=(4d0*ex(i-1,j,k)-ex(i-2,j,k))/3d0
      elseif(bcMx.eq.-2) then
        !apply open, Horiuchi01
        ey(i,j,k)=(5d0*ey(i-1,j,k)-4d0*ey(i-2,j,k)+ey(i-3,j,k))/2d0
        ez(i,j,k)=(4d0*ez(i-1,j,k)-ez(i-2,j,k))/3d0
        ex(i,j,k)=(4d0*ex(i-1,j,k)-ex(i-2,j,k))/3d0
      else
        !apply open ... Pritchett ..
        ey(i,j,k)=(4d0*ey(i-1,j,k)-ey(i-2,j,k))/3d0
        ez(i,j,k)=0d0
        ex(i,j,k)= (ex0(i,j,k)-susxy(i,j,k)*ey(i,j,k)-susxz(i,j,k)&
               *ez(i,j,k) -jxtilde(i,j,k)*cntr*dt*1.d0)/susxx(i,j,k)*1.d0
	  end if

	  end if
	  enddo
	  enddo
!
!	Boundaries in Y
!

	j=2
  do i=2,nxp
      do k=2,nzp
         if (periodicy.or.bcMy.eq.-1) then
            ! do nothing
         else

            if(bcMy.eq.1) then
!   		   ! apply conductor
               ex(i,j,k)=ex_ext
               ey(i,j,k)=0d0 !needs to be corrected
               ez(i,j,k)=ez_ext
            else
               !apply mirror
               ex(i,j,k)=(4d0*ex(i,j+1,k)-ex(i,j+2,k))/3d0
               ey(i,j,k)=0d0
               ez(i,j,k)=(4d0*ez(i,j+1,k)-ez(i,j+2,k))/3d0
            end if

         end if
      enddo
   enddo


	j=nyp

 do i=2,nxp
      do k=2,nzp
         if (periodicy.or.bcMy.eq.-1) then

	           ex(i,j,k)=ex(i,2,k)
  	           ey(i,j,k)=ey(i,2,k)
  	           ez(i,j,k)=ez(i,2,k)

         else

            if(bcMy.eq.1) then
!              ! apply conductor
               ex(i,j,k)=ex_ext
               ey(i,j,k)=0d0 !needs to be corrected
               ez(i,j,k)=ez_ext
            else
               !apply mirror
               ex(i,j,k)=(4d0*ex(i,j-1,k)-ex(i,j-2,k))/3d0
               ey(i,j,k)=0d0
               ez(i,j,k)=(4d0*ez(i,j-1,k)-ez(i,j-2,k))/3d0
            end if

         end if
      enddo
   enddo
!
!	Boundaries in Z
!
!	min boundary in z

   k=2
   do i=2,nxp
      do j=2,nyp
         if (periodicz) then
            ! do nothing
         else

            if(bcMz.eq.1) then
               ! apply conductor
               bavg2=bx(i,j,k)**2+by(i,j,k)**2+bz(i,j,k)**2
               g_x=sin(pi*x(i,j,k)/dx/dble(ibar))**2
               if(a_m.gt.bx_ini.or.a_m.lt.bx_ini*0.000001)then
                  v_int=2.d0*a_m*o_m*tanh(o_m*t)/cosh(o_m*t)**2
                  !meaning Ey_ext is big or zero, that is classical Newton or nothing
               else
                  v_int=a_m*tanh(o_m*t*5d0)**2
                  !for the constant Ey case, Ey=0.1, 0.05 or 0.01 of V_A=bx_ini
               endif
               if(shock) then
                  g_x=1.d0
                  v_int=ey_ext
               end if
!
               ex(i,j,k)=ex_ext
               ey(i,j,k)=(v_int*g_x + 0.3849d0*2.d0*a_m*o_m*tanh(o_m*t*1d0)**2*0d0)&
                    *bavg2/(abs(bx(i,j,k))+1.d-10)
               !ey(i,j,k)=v_int*g_x*(abs(bx(i,j,k))+1.d-10)
               ez(i,j,k) = (ez0(i,j,k)-suszx(i,j,k)*ex(i,j,k)-suszy(i,j,k)&
                    *ey(i,j,k) -jztilde(i,j,k)*cntr*dt*1.d0)/suszz(i,j,k)*1.d0
            elseif(bcMz.eq.2) then
               !apply mirror
               ex(i,j,k)=(4d0*ex(i,j,k+1)-ex(i,j,k+2))/3d0
               ey(i,j,k)=(4d0*ey(i,j,k+1)-ey(i,j,k+2))/3d0
               ez(i,j,k)=0d0
            elseif(bcMz.eq.0) then
               !apply open
               ex(i,j,k)=(4d0*ex(i,j,k+1)-ex(i,j,k+2))/3d0
               ey(i,j,k)=(4d0*ey(i,j,k+1)-ey(i,j,k+2))/3d0
               ez(i,j,k) = (ez0(i,j,k)-suszx(i,j,k)*ex(i,j,k)-suszy(i,j,k)&
                    *ey(i,j,k) -jztilde(i,j,k)*cntr*dt*1.d0)/suszz(i,j,k)*1.d0
            elseif(bcMz.eq.-1) then
               !apply open
               ex(i,j,k)=(4d0*ex(i,j,k+1)-ex(i,j,k+2))/3d0
               ey(i,j,k)=(4d0*ey(i,j,k+1)-ey(i,j,k+2))/3d0
               ez(i,j,k) = (ez0(i,j,k)-suszx(i,j,k)*ex(i,j,k)-suszy(i,j,k)&
                    *ey(i,j,k) -jztilde(i,j,k)*cntr*dt*0.d0)/suszz(i,j,k)*1.d0
            else
               !apply open ... Pritchett..
               ex(i,j,k)=0d0
               ey(i,j,k)=(4d0*ey(i,j,k+1)-ey(i,j,k+2))/3d0
               ez(i,j,k)= (ez0(i,j,k)-suszx(i,j,k)*ex(i,j,k)-suszy(i,j,k)&
                    *ey(i,j,k) -jztilde(i,j,k)*cntr*dt*1.d0)/suszz(i,j,k)*1.d0
            end if

         end if
      enddo
   enddo
!
!	end

	k=nzp
      do j=2,nyp
	  do i=2,nxp
	  if (periodicz) then
	  ex(i,j,k)=ex(i,j,2)
  	  ey(i,j,k)=ey(i,j,2)
  	  ez(i,j,k)=ez(i,j,2)
	  else

	  if(bcMz.eq.1) then
	    ! apply conductor
	    bavg2=bx(i,j,k)**2+by(i,j,k)**2+bz(i,j,k)**2
        g_x=sin(pi*x(i,j,k)/dx/dble(ibar))**2
        if(a_m.gt.bx_ini.or.a_m.lt.bx_ini*0.000001)then
           v_int=2.d0*a_m*o_m*tanh(o_m*t)/cosh(o_m*t)**2
! meaning Ey_ext is big or zero, that is classical Newton or nothing
        else
           v_int=a_m*tanh(o_m*t*5d0)**2
!! for the constant Ey case, Ey=0.1, 0.05 or 0.01 of V_A=bx_ini
        endif
        if(shock) then
           g_x=1.d0
	       v_int=ey_ext
        end if

	    ex(i,j,k)=ex_ext
	    ey(i,j,k)=(v_int*g_x + 0.3849d0*2.d0*a_m*o_m*tanh(o_m*t*1d0)**2*0d0)&
                 *bavg2/(abs(bx(i,j,k))+1.d-10)
        !ey(i,j,k)=v_int*g_x*(abs(bx(i,j,k))+1.d-10)
        ez(i,j,k) = (ez0(i,j,k)-suszx(i,j,k)*ex(i,j,k)-suszy(i,j,k)&
               *ey(i,j,k) -jztilde(i,j,k)*cntr*dt*1.d0)/suszz(i,j,k)*1.d0
	  elseif (bcMz.eq.2) then
	    !apply mirror
        ex(i,j,k)=(4d0*ex(i,j,k-1)-ex(i,j,k-2))/3d0
        ey(i,j,k)=(4d0*ey(i,j,k-1)-ey(i,j,k-2))/3d0
   	    ez(i,j,k)=0d0
      elseif(bcMz.eq.0) then
        !apply open
        ex(i,j,k)=(4d0*ex(i,j,k-1)-ex(i,j,k-2))/3d0
        ey(i,j,k)=(4d0*ey(i,j,k-1)-ey(i,j,k-2))/3d0
        ez(i,j,k) = (ez0(i,j,k)-suszx(i,j,k)*ex(i,j,k)-suszy(i,j,k)&
               *ey(i,j,k) -jztilde(i,j,k)*cntr*dt*1.d0)/suszz(i,j,k)*1.d0
      elseif(bcMz.eq.-1) then
        !apply open
        ex(i,j,k)=(4d0*ex(i,j,k-1)-ex(i,j,k-2))/3d0
        ey(i,j,k)=(4d0*ey(i,j,k-1)-ey(i,j,k-2))/3d0
        ez(i,j,k) = (ez0(i,j,k)-suszx(i,j,k)*ex(i,j,k)-suszy(i,j,k)&
               *ey(i,j,k) -jztilde(i,j,k)*cntr*dt*0.d0)/suszz(i,j,k)*1.d0
      else
        !apply open ... Pritchett ..
        ex(i,j,k)=0d0
        ey(i,j,k)=(4d0*ey(i,j,k-1)-ey(i,j,k-2))/3d0
        ez(i,j,k)= (ez0(i,j,k)-suszx(i,j,k)*ex(i,j,k)-suszy(i,j,k)&
             *ey(i,j,k) -jztilde(i,j,k)*cntr*dt*1.d0)/suszz(i,j,k)*1.d0
	  end if

	  end if
	  enddo
	  enddo

      end subroutine bc_maxwell








